LIGO-P060063-02-Z 



Binary system delays and timing noise in searches for gravitational waves from known 

pulsars 

Matthew Pitkin and Graham Woan 1 

1 Department of Physics and Astronomy, University of Glasgow, University Avenue, Glasgow, G12 8QQ, UI^ 

(Dated: February 7, 2008) 

The majority of fast millisecond pulsars are in binary systems, so that any periodic signal they 
emit is modulated by both Doppler and relativistic effects. Here we show how well-established 
binary models can be used to account for these effects in searches for gravitational waves from 
known pulsars within binary systems. A seperate issue affecting certain pulsar signals is that of 
timing noise and we show how this, with particular reference to the Crab pulsar, can be compensated 
for by using regularly updated timing ephemerides. 
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I. INTRODUCTION 

Neutron stars are thought to be strong candidates 
for the emission of detectable continuous gravitational 
waves including the 1627 pulsars currently discov- 
ered 1 . The majority of these pulsars have been discov- 
ered through dedicated radio surveys. Surveys are on- 
going, but estimates of the number of active pulsars in 
the galaxy can be made by inference from the current 
population, taking into account biasing from selection 
effects, and the supernova rate. Estimates give values of 
~ 200 000 active pulsars within our galaxy (see Ref. [H). 

Pulsars are found in a wide range of environments. 
Some are directly associated with the supernova rem- 
nants (SNRs) in which they were born. These are typ- 
ically young pulsars whose birth velocity has not yet 
caused a large displacement from the remnant and the 
SNR has not dissipated into the interstellar medium 
(ISM). Some pulsars are found in binary systems as the 
companions of a wide range of astronomical bodies from 
planets, through main sequence stars, to white dwarfs 
and other neutron stars. The fastest, 'millisecond', pul- 
sars (pulsars with rotation periods of < 10 ms) are usu- 
ally found within binary systems, and often within globu- 
lar clusters, their rapid rotation rate a consequence of be- 
ing spun-up by accretion of material from a stellar com- 
panion Q . Pulsars are also seen without any association, 
and in this paper we will classify any pulsar not in a bi- 
nary system as isolated. 

Pulsars are generally seen to spin-down as they lose 
rotational energy through a variety of emission mecha- 
nisms, but the primary loss mechanism is thought to be 
magnetic dipole radiation. Other potential mechanisms 
include particle acceleration and gravitational radiation. 
Whatever the mechanisms at work, the rotational phase 
evolution of a pulsar can generally be well described by 



a short Taylor expansion, 



(f>(T) = O + 27r<j ^ (r - *o) + -MT-to) 2 



+ gi>b(T-to) 3 + 



(1.1) 



where 4>o is the initial phase, v$ and its time derivatives 
are the pulsar frequency and spin-down coefficients at an 
epoch to, and T is the time in a frame comoving with the 
pulsar. For the vast majority of pulsars the value of v is 
very small and v is unmcasurable or swamped by timing 
noise (see fllVj) . Note that gravitational waves emitted 
from a triaxial, non-precessing, neutron star come from 
the quadrupolar component of the rotating body and so 
will have exactly twice the phase evolution described by 
Eq.O 

The expected gravitational wave signal from a triaxial 
neutron star is given by [|| 



h(t) = -F + (t;tp)h (l + cos 2 t) cos 2^) 
+F X (t; ip)ho cos l sin 2<p(t), 



(1.2) 



where cf>(t) is that given in Eq. (|1.1|) . F + and F x arc the 
detector beam patterns for the plus and cross polarisa- 
tions of the gravitational waves, ip is the wave polarisa- 
tion angle, and i is the angle between the rotation axis of 
the pulsar and the line-of-sight. For a gravitational wave 
signal impinging on the Earth the signal arrival time at 
the detector, t, will be modulated by Doppler, time de- 
lay and relativistic effects caused by the motions of the 
Earth and other bodies in the solar system, so we want 
a reference frame which will not be affected by these. In 
general, for isolated pulsars 2 , such a frame is the solar 
system barycentre (SSB). To convert from t to the time 
at the SSB, tb, we must include a series of time correc- 
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1 As given by the Australia Telescope National Facility - ATNF 
- online pulsar catalogue 0, Q| as of 6 th Nov 2006, from which 
time all subsequent pulsar numbers will be taken. 
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tions, 

t b = t + St = t+ — + A E@ +A Se , (1.3) 

where r is the position of the detector with respect to the 
solar system barycentre (SSB), n is the unit vector point- 
ing to the pulsar, A^ is the special rclativistic Einstein 
delay, and Ag is the general relativistic Shapiro delay 
(see Ref. Q for definitions of these delay terms). For 
isolated pulsars we can therefore set T = t b in Eq. (|1 . ip . 
For pulsars in binary systems there will be additional 
time delays as discussed in §1111 

In this paper we describe the addition of such binary 
system time delays with respect to current searches for 
gravitational waves from known pulsars. In a seperate 
issue we also discuss how pulsar timing noise can cause 
deviations from the simple phase model in Eq. 11.11 and 
how to account for this in the analysis. 



II. SEARCH METHOD 

Searches for gravitational waves from a selection of 
known pulsars have been performed using data from the 
LIGO (lo| and GEO 600 [ll[ gravitational wave detec- 
tors from the four science runs (Sl-4), which have taken 
place since late 2002 [H [H Q. The search method 
used is outlined briefly here, but is described more fully 
in Refs. HQ. 

All current data from these detectors is sampled at 
16 384 Hz, giving a range of 0-8192 Hz available for 
searches, although in practice this is limited to the most 
sensitive region of the detectors between « 50-4000 Hz. 
The rotation frequencies of known pulsars are, of course, 
known very precisely from radio and/or X-ray observa- 
tions, so the vast majority of this frequency space is re- 
dundant and the speed of any search can be increased by 
removing it. Knowledge of the pulsar parameters allows 
us to perform a complex heterodyne on the data, using 
the precise pulsar phase evolution, and down-sample it 
to Hz. Following this heterodyne a Bayesian parame- 
ter estimation for the unknown pulsar parameters /io, ip, 
cos l and 4>o can be performed on the massively reduced 
data set. 

The method relies on the accuracy of the signal model, 
and it is very important that the phase evolution in the 
heterodyne is sufficiently good. Any drift from the true 
pulsar signal phase could nullify the search if it becomes 
too severe. This paper will discuss how to ensure that 
this phase model is sufficiently accurate for pulsars in bi- 
nary systems and those strongly affected by timing noise. 
Most recently, these methods have been used in Ref. [l4| 
to obtain limits on the gravitational wave emission from 
78 pulsars using data from the LIGO and GEO 600 S3 
and S4 runs. 



III. PULSARS IN BINARY SYSTEMS 

Of the 1627 pulsars in the ATNF catalogue, 124 are in 
binary systems. The first of these to be discovered was 
the highly relativistic pulsar PSR J1915+1606 found by 
Hulse and Taylor in 1974 [lj]. Of these 124 pulsars, 98 
have spin frequencies greater than 25 Hz (out of a total 
163 isolated and binary pulsars with frequencies greater 
than 25 Hz), and therefore gravitational wave frequen- 
cies > 50 Hz, putting them into the sensitive band of the 
LIGO detectors. Indeed, for the reason stated in jjj the 
majority of millisecond pulsars are in binary systems. 

A. Pulsar timing 

A brief discussion on how pulsar timing information is 
obtained is relevant here. The majority of pulsars have 
been both discovered and monitored in radio. Pulsar 
surveys, discussed in more detail in Ref. [l6j ]. use Fourier 
transform methods to look for periodic signals in the ra- 
dio data, taking into account the effect of interstellar dis- 
persion across the receiver band. Once the pulsar period 
has been determined, the radio time series data can be 
folded with this cadence to build up the signal-to-noise 
ratio of a mean pulse. Once a stable pulse is obtained 3 
the time of arrival (TOA) can be measured at the peak of 
the pulse. These pulse times can then be used to extract 
more precise information about the pulsar parameters, 
including its position and frequency parameters. 

The most prevalent tool used to fit timing measure- 
ments is the Tempo software suite [l7[, however others 
have been developed, including psrtime [l8[ at Jodrell 
Bank Observatory. Tempo requires precise solar system 
ephemerides, containing the positions and velocities of 
the major solar system bodies, to convert TOAs at a de- 
tector to the rest frame of the pulsar. It computes the 
pulsar phase at each TOA, </>(Ti), over the range of pulsar 
parameters (a, 5, v, z>, etc), and uses a x 2 goodness-of-fit 
statistic to determine the best model via minimization. 
A starting point for the fit is obtained through a rough 
knowledge of the position and frequency from the initial 
discovery, but it can still be quite complex as there can 
be many other parameters that could be contributing. 
Below it is seen how a pulsar in a binary system requires 
a complex model with many more parameters than an 
isolated object. 

B. Binary pulsar timing 

The majority of pulsars within the ground-based grav- 
itational wave frequency band are in binary systems, so 



Individual pulses can vary in shape, but the summation of many 
gives a generally stable pulse shape. 
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we must consider the effects of binary timing carefully. 
Techniques to create filters to match a binary signal (in 
the frequency domain) have been described before, for 
example in Ref. [l9l ]. and these have been applied to 
the search f or g ravitational waves from the neutron star 
in Sco X-l [2(|, but there has been no equivalent treat- 
ment in the time domain, suitable for the search outlined 
above. Eq. (|1.3|) shows the timing corrections needed to 
take account of Doppler and relativistic delays of a signal 
and transform it to the SSB. Any constant Doppler de- 
lays from the pulsar's actual motion relative to the SSB 
are unimportant, and the SSB frame can be considered 
as the rest frame of the pulsar. For a pulsar in a bi- 
nary system however, its motion within the system will 
need to be taken into account with a transform from the 
binary system barycentre to the pulsar proper time. 

The basic transformation and binary models below are 
summarised in Ref. 0] and used in the pulsar timing 
program Tempo [l7| ■ The transformation from SSB time 
if, to pulsar proper time T follows the form of Eq. (|1.3p 
and is 

t 6 -i =T + A R + A E + A s +A A , (3.1) 

where Apj is the Roemer time delay giving the propa- 
gation time across the binary orbit; Ag is the Einstein 
delay and gives gravitational rcdshift and time dilation 
corrections; Ag is the Shapiro delay which gives the grav- 
itational propagation delay due to the signal propagation 
through the curved space-time of the companion; and Aa 
is the aberration delay caused by the pulsar's rotation. 

The majority of binary orbits can be well-modelled as 
Keplcrian. Keplcrian orbits are defined by five param- 
eters, To - the time of periastron (closest approach in 
the binary orbit); to - the longitude of periastron; F\ 
- the orbital period; e - the orbital eccentricity (where 
e = 1^/(1 — b 2 / a 2 ) and a and b are the semi- major and 
semi- minor axis of the orbital ellipse respectively); and 
x = (asini)/c is the projected semi-major axis, with i 
being the orbital inclination. For some of the most ex- 
treme binary systems, with rapidly-spinning pulsars in 
tight orbits, the maximum (gravitational wave) orbital 
Doppler frequency shift can be up to ~ 0.1 Hz. For ex- 
ample, PSR J0737-3039A with x = 1.415 light sec, P b = 
0.10225 days and u gw = 88.11Hz gives Ai/ gw = 0.089 Hz. 

We will now consider the three most commonly used 
models used to characterise the TOA of pulses from pul- 
sars in binary systems. 

1. Blandford-Teukolsky model 

The Blandford and Tcukolsky model [H| (BT) makes 
no assumptions about the correct theory of gravity. In- 
stead, it assumes a simple Keplerian orbit with slow 
precession, into which additional relativistic effects have 
been added. Other phenomena can be taken into ac- 
count through time derivatives of the four main orbital 
elements, excluding Tq. The BT model has been used to 



fit data for 53 of the binary pulsars with v > 25 Hz, and 
is the most common model used. Exceptionally, one of 
these systems is modelled using the Tempo model BT2P 
which accommodates three orbits, the first of which can 
be relativistic, but the second and third are Keplerian. 
The system is a multiple system, described in Ref. [22j], 
in which three, or possibly four, planets orbit the pulsar. 
Although these additional orbits complicate the above 
equations, the standard BT model is sufficient for our 
purposes. 



2. Low eccentricity model 

The second most common model used in fitting ra- 
dio observations of binaries is the low eccentricity model 
(called ELL1 in Tempo) developed in Ref. [lU- It is used 
as a fit for pulsars in very low eccentricity orbits (e ~ 0) 
and is the appropriate model for 38 of our target pulsars 
with v > 25 Hz. An almost circular orbit makes To an d oj 
nearly degenerate, so these parameters, along with e, are 
replaced with the non-covariant parameters of the time 
of the ascending node of the orbit (T asc = T — u)P b /2ir) 
and the first and second Laplace-Lagrange parameters 
?y = esinw and k = ecosw. The time delays for this 
model are defined in [23[ and Tempo. 



3. D amour- Deruelle model 

The third most common model is that of Damour and 
Deruelle (DD) [24|. This model uses a method for solv- 
ing the relativistic two-body problem to post-Newtonian 
order and is valid under very general assumptions about 
the nature of gravity in strong field regimes. It is useful 
for highly relativistic systems, although in only mildly 
relativistic systems this model should be no different to 
the BT model. There are seven pulsars with v > 25 Hz 
in the ATNF catalogue that use this model. This model 
is again summarised in [3]. 



C. Comparison with TEMPO 

The above three models are all implemented in the pul- 
sar timing software package Tempo. In our search for 
gravitational waves from binary systems we also require 
these additional time corrections to correctly calculate 
the phase of the pulsar for heterodyning. Code to cal- 
culate the binary time delays for each model has been 
adapted from the Tempo counterparts and is available 
under CVS in the LIGO Algorithm Library (LAL) repos- 
itory [25[ . Some consistency tests have been performed 
between the two codes, which are described below. It 
is acknowledged that Tempo has some uncertainties in 
its timing models due to simplifying assumptions, and 
these limit its accuracy to ^100 ns. These effects are dis- 
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cussed in Refs. (2(| [27}, but as errors at this level have 
no significant impact on our search we will neglect them. 

1. PSR J1012+5307 

We performed a validation check on the new LAL code 
by demodulating radio data from a known pulsar. With 
no known gravitational calibrators, this is a crucial step 
in the development of any gravitational wave detection 
code. A set of TOAs for PSR J1012+5307 obtained with 
the Effelsberg 100 m radio telescope in Bonn, Germany, 
was kindly supplied by Michael Kramer for this purpose. 
This pulsar has the second most circular orbit known 
and is therefore fitted by the ELL1 model. The data 
intermittently spanned just over 5 years from 2 nd Jan- 
uary 2000 to 12 th February 2005 and comprised TOAs 
(in Modified Julian Date format, where MJD = Julian 
Date-2400000.5), together with the TEMPO-derived pul- 
sar parameters (Table [XJ) . Some minor transformations 

TABLE I: The parameters of PSR J1012+5307. Values are 
quoted with la errors on the final digit in brackets. 



PSR. .11012+5307 



a 


10 h 12 m 33 a . 43368(1) 


5 


53°07'02". 5880(2) 


PMRA 


2.38(3) mas/yr 


PMDEC 


— 25.35(5) mas/yr 


V 


190.267837621884(9) Hz 


if 


-6.2022(2) X 10" 16 Hz/s 


V 


2.0(3) x 10~ 27 Hz/s 2 


Frequency epoch 


MJD 50700 


Dispersion measure 


9.0233(7) cm" 3 pc 


Observing Frequency 


1408.6 MHz 


Binary model 


ELL1 


X 


0.581817(1) s 


p b 


0.6046727136(2) days 


T aac 


MJD 50700.0816289(4) 


V 


7(4) x 10~ 7 


K 


-1(40) x 10~ 8 



are necessary to convert TOAs measured at the telescope 
to the GPS time stamps used in gravitational wave data 
analysis software. First the raw TOAs are corrected for 
the drifts between the hydrogen maser clock at Effelsberg 
and coordinated Universal Time of the National Institute 
of Science and Technology UTC(NIST) reference. This 
correction (supplied with the data) was typically a few 
microseconds. The difference between UTC(NIST) and 
UTC has been less than ±100 ns since 6 th July, 1994 (H 
and was neglected for this work. The conversion between 
the time scales therefore becomes 

icps = (iuTC(MjD) - 44244 days) x 86400 s + L, (3.2) 

where the 44244 corresponds to the MJD of the GPS 
time epoch (1 st January, 1980) and L is the accumulated 
number of leap seconds included in the definition of UTC. 
For the time-span of these TOAs, L = 13. 



The TOAs can now be corrected for interstellar dis- 
persion time delay (l6l | 

A< disp = 4.149 x 10 3 MHz 2 pc" 1 cm 3 s x DM// 2 s, (3.3) 

where DM is the dispersion measure in cm -3 pc and / is 
the radio observation frequency in MHz (sec Table Q] for 
values). This correction is subtracted from the TOAs to 
give observations at infinite frequency with no dispersion. 

One of the major differences between our binary time 
domain code and the Tempo code is the time system 
used. All epochs in Tempo arc defined as MJD Barycen- 
tric Dynamical Time (TDB - a timescale generally used 
for ephemerides referenced to the solar solar barycentre) 
whereas the general reference time for our gravitational 
wave data is GPS time. Epochs therefore have to be con- 
verted to GPS time on the TDB timescale. This TDB 
timescale is related to Terrestrial Time (TT - formerly 
Terrestrial Dynamical Time TDT), which represents a 
time consistent with relativity for an observer on the 
Earth's surface, by a small factor, TDB = TT + St, no 
greater than a couple of milliseconds and given by 

St = 0.001 658 s x sin $ + 0.000 014 s x sin 2$, (3.4) 

where $ = 357.53° + 0.985 600 28°(MJD - 51544.5) is 
the mean anomaly, or phase, of the Earth's orbit at the 
given Modified Julian Date [29]. TT is offset from In- 
ternational Atomic Time (TAI), so that TT = TAI + 
32.184 seconds 4 . The conversion thus goes ^tdb(gps) = 
(^tdb(mjd) - 44244) - 51.184 - St, where the 51.184s 
comprises the 32.184 s difference between TT and TAI 
and 19 second difference between TAI and GPS. 

The LAL code to calculate the SSB time delay uses 
the pulsar's position, the telescope position and a solar 
system ephemeris [301 an d was use d for each pulsar TOA 
to correct to the time at the SSB. This code was writ- 
ten by Curt Cutler and has been independently tested 
against Tempo [1, [l2j showing no more than 4 //s differ- 
ence between the two. For a pulsar with a gravitational 
wave frequency at 1 kHz a 4 fis timing error would give 
a phase offset between any signal and our model tem- 
plate of 0.025 rads. This would reduce our sensitivity to 
a signal by of order 1 — cos 0.025 « 3x 10 -4 , which is 
negligible. 

Once corrected to the SSB the TOAs are further cor- 
rected to the pulsar proper time by calculating the time 
delays in the binary system using the binary system pa- 
rameters (see Table Q}. The binary and solar system time 
delays for a selection of TOAs covering part of the binary 
orbit are shown in Fig. [1] 

Once these corrections to the TOAs are applied we can 
compare the LAL barycentring codes with Tempo by in- 
serting the TEMPO-derived pulsar parameters into the 



4 There are many definitions of time used in astronomy and very 
careful attention of which one is being used and how to convert 
between them is essential when high precision timings are being 
made. A good guide for these is |29|| . 
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-0.5 - 
6.32C 

3 65 . 4| — 
365.2 - 
365 - 
364.8- 
364.6- 



1 6.3201 6.3201 6.3202 6.3202 6.3203 6.320 
GPS time (sees) 



04 6.3205 

x 10 B 



01 6.3201 6.3201 6.3202 6.3202 6.3203 6.320 
GPS time (sees) 



04 6.3205 
x 10 s 



first TOA (years) 



FIG. 3: The modulus of the pulsar phase at each TO A over 
a 5 year period corrected for the binary delay. 



FIG. 1: The binary and solar system time delays calculated 
for PSR J1012+5307 over a part of a binary orbit. 



LAL baryecntring routines and examining the predicted 
phase at the corrected TOAs. TOAs converted incor- 
rectly by the baryecntring codes would show up as a 
phase drift. The phase at each TOA was calculated 
using the supplied frequency and frequency derivatives 
in Eq. (jl.lj) . with 4>q = and the frequency epoch as 
to. Fig. [5] demonstrates the serious effect of neglecting 
these binary time delays. In contrast, Fig. [3] shows how 




12 3 4 

Time from first TOA (years) 



FIG. 2: The modulus of the pulsar phase at each TOA over 
a 5 year period with no binary time delays included. 



the TOAs barycentred using our code stay well in phase 
over the observation time when the binary delays are in- 
cluded, with a residual phase slope of ~ 0.04rads/yr. A 
yearly periodicity is also present possibly showing up the 
slight difference in the LAL solar system barycentring 
code and Tempo, although these effects are at a very 
low level. Several points clearly show large phase resid- 
uals and correspond to times when the level of noise on 
the TOA measurements was high. 

The parameters for PSR J1012+5307 were generated 
using the ELL1 model, so the above test only checked 
the ELL1 code. We can also check the two other models 
by converting T asc to T and the Laplacc-Lagrange pa- 
rameters k and r\ to e = \J (k 2 + rj 2 ) and lu. This pulsar 



has a low eccentricity, so T can be set equal to T asc and 
e and oj set to zero for practical purposes. Doing this 
we can again produce the phase plots for the BT and 
DD models (Fig. [4|. The phase is again well described 



Time from first TOA (years) 



Time from first TOA (years) 



FIG. 4: The modulus of the pulsar phase at each TOA over 
a 5 year period for the BT and DD models, corrected for the 
binary delay. 



for these two models, with the slope and periodicity still 
present. This suggests that the slope and periodicities 
are not caused by the binary timing correction code (as 
each model is independent), but may be a results of slight 
errors in the other timing corrections, the solar system 
barycentring code, the ephemerides, or the parameters 
used. 

The effects of inaccuracies in these binary parameters 
can be shown by offsetting one from its true value. In 
Fig-E]the true value of the T asc parameter has been offset 
by 5, 10 and 20 seconds and the phase at each TOA 
recalculated. Even for a 5 s mismatch we start to depart 
from the true phase by up to 0.5 radians. This would 
not be disastrous for the analysis but would degrade the 
search. For a 20 s offset we start to see phase errors 
of up to 2 radians. At this level, all sensitivity is lost. 
In real analyses as in Ref. [l4| the effect of any errors 
on the parameters, as to whether they could cause our 
heterodyne phase to depart significantly from the signal 
phase, is thoroughly checked. 



5 



C Comparison with TEMPO 



IV THE PROBLEM OF TIMING NOISE 



20 sees 
10 sees 
5 sees 



hk:.. ... .■ ..; ... • ii V. j 



v i , 



12 3' 
Time from first TOA (years) 



FIG. 5: The modulus of the pulsar phase at each TOA over 
a 5 year period with T asc offset from its true value. 



2. Direct check against TEMPO 

In common with the solar system baryecntring code 
@j EH j the binary timing code was tested directly against 
Tempo. Tempo can be run in predictive mode, to use a 
set of pulsar parameters to predict the pulsar phase over 
a period of time. This predicted phase can be then be 
compared with that calculated using our binary timing 
code. This was done for each model with a set of 100 
randomly generated binary pulsar systems over a period 
of 100 days. The detector location was set to be at the 
SSB, so the solar system time delay errors would not 
be included. Histograms of the time residuals between 
the codes are shown in Fig. [6] for each model. These 
show that the time difference between the two codes is 
generally less than ±1 fj,s, which is sufficiently small to 
ensure any signal and template remain in phase over the 
frequency range considered. 



IV. THE PROBLEM OF TIMING NOISE 

A rather different issue that if unaddressed could po- 
tentially cause problems for our analysis is that of timing 
noise. 

Pulsars are generally very stable over periods of sev- 
eral days, but there are phenomena which can cause de- 
viations in this timing stability. With the very high ac- 
curacy of pulsar timing any random timing irregularities 
will start to become evident. These phenomena show 
up as glitches and timing noise. Timing noise has been 
known about since the early days of pulsar observations 
and represents a random walk in phase, frequency or fre- 
quency derivative of the pulsar about the regular spin- 
down model given in Eq. (jTTTJ) [3l[ . If our search method 
cannot take account of these phase deviations they could 
reduce our sensitivity or even invalidate the search. It is 
therefore important to see firstly if timing noise produces 
large enough phase deviations to cause a problem for our 



method and if so how it can be countered. 

The strength of this effect has been quantitatively de- 
fined in Ref. [3l[ as the activity parameter A, as refer- 
enced to that of the Crab pulsar, and in Ref. [32[ as the 
stability parameter As. There is however no real con- 
sensus on how to quantitatively define a measure of the 
level of timing noise, with the magnitude and sign of P 
maybe providing the measure which includes the least 
other assumptions. A thorough study of timing noise, 
comparing and contrasting the various measures used, 
is given in Ref. [H] (also see Refs. [H, [35[). There is 
a definite correlation between these parameters (A and 
Ag) and the pulsar's spin-down rate, therefore possibly 
the pulsar's age. Young pulsars, like the Crab pulsar, 
generally show the most timing noise activity. 

The Crab pulsar is the youngest known pulsar targeted 
by current gravitational wave detectors and it is impor- 
tant to considered how timing noise may be countered 
for this pulsar. A method was first proposed in Ref. 13611 
and has been used in the analysis of Abbott et al. [13], 
but has not previously been described in detail. A brief 
look into the effects of timing noise for other pulsars will 
be discussed later. 



A. Timing noise in the Crab pulsar 

The pulsar (J0534+2200) in the Crab nebula (Ml) has 
undergone intense study since its discovery in 1968. Its 
parameters are given in Table [TTJ It has an observed 



TABLE II: The parameters of the Crab pulsar as calculated 
from the Jodrell Bank monthly ephemeris [39l ] . 



PSRJ0534+2200 



Right ascension a 


05 h 34 m 31 s .973 


Declination 5 


22°00'52".06 


proper motion in a 


— 13 mas/yr 


proper motion in 5 


7 mas/yr 


Position epoch 


MJD 40675 


V 


29.7670971390 Hz 


V 


-3.72633xlO~ 10 Hz/s 


V 


1.1429 x!0~ 20 Hz/s 2 


Frequency epoch 


MJD 53993 


Distance 


2.0 kpc 



age of 951 years (the formation of the Crab nebula is 
associated with a supernova observed in AD 1054) and a 
spin-down age of —v/2v = P/2P = 1250years. Analyses 
of long-term timing observation of the Crab pulsar are 
given in Refs. [13, HH- These analyses show some of 
the timing features which make the Crab pulsar such an 
interesting object: the timing noise and glitches. 

Since 1982 there has been a regular monitoring pro- 
gram of the Crab pulsar at Jodrell Bank Observatory, 
and ti ming ephemerides from this are publicly available 
online [39(. The ephemeris gives the pulsar frequency 
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timing residual (|rs) timing resdiual (|Is) timing residual (JIs) 

FIG. 6: Timing residuals between the pulsar phase as predicted by Tempo and that computed with our binary code for 100 
random pulsars for each binary model. 



and frequency derivative and associated errors, and the 
associated epoch. The epochs, generally given on the 
15 th of each month, represent the time of the peak of the 
first pulse after midnight on that day. They therefore 
represent zero of modulus phase of the electromagnetic 
pulse. Notes are given in the event of a timing irregular- 
ity or glitch being observed. The Crab's timing noise can 
be clearly seen in the ephemeris once a best-fit quadratic 
timing model (Eq. has been subtracted (see Fig. [7]). 

The section of data used was chosen to be free of glitches 
as these can be much larger than any timing noise fre- 
quency deviations. Fig. [7] compares well with that given 




45000 45500 46000 

Modified Julian Date (days) 



FIG. 7: The timing noise in the frequency of the Crab pulsar 
after removing a quadratic fit to the frequency as given in the 
Jodrell Bank ephemeris. 



in Ref. [37J, although some difference can be expected 
due to the different lengths of data and epochs used in 
the fitting. It can be seen that on scales of several months 
there is quite a large variation in the timing residual. It 



is shown in Ref. [37| that on smaller time scales the varia- 
tion is far smoother. It is important to consider whether 
timing noise at this level will cause a loss in sensitivity in 
any search for gravitational wave emission from the Crab 
pulsar. Jones [4(| constructs a dccoherence timescale, 
Tdccohcrencej defined as the time over which the timing 
noise will cause the phase to deviate by 1 radian from 
the second order Taylor expansion of phase. This makes 
use of the "activity parameter" and is calculated to be 
~ 2.6 yr for the Crab pulsar. Though useful as an indi- 
cator, this statistic does not take account of permanent 
changes to spin-down caused by glitches or other secular 
variations in the pulsar. 

Although it is reasonable that a third order fit to the 
entire data would improve the model, there is no need. 
The Crab pulsar ephemeris provides timing every month, 
which is sufficient to track the phase excursions in the 
timing noise. All that is necessary is that these be in- 
terpolated between the ephemeris times. By using the 
phase, frequency and frequency derivative for each entry 
in the ephemeris as boundary conditions to a set of simul- 
taneous equations the full phase evolution between each 
month can be calculated, giving a fifth order polynomial, 

foth(r) = 0o + 2tt jV (T - t ) + h) Q (T - t ) 2 + 

UhiT - t ) 3 + ^ (T - t ) 4 + ^o(T - to) 5 |(4.1) 

Indeed, for much of the time even this method is un- 
necessarily complicated and a simple linear interpolation 
between months would be sufficient. 

These corrections can be included in the method de- 
tailed in iJTT] as an extra heterodyne step [36| . The initial 
heterodyne (described in fJTTJ) , uses a third order fit to 
the the phase with values of v and i> taken from the 
ephemeris at the closest time before the timestamp on 
the data to be analysed, and v taken from the ATNF 
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pulsar catalogue value given in Table IIIII Then, assum- 
ing that any gravitational wave signal would show the 
same timing noise (see Refs. [1(3] and [3] for discussions 
of this), we apply a second heterodyne to the data 
using the phase difference between Eqs. and (|4.ip 



(T)] 



(4.2) 



where the factor of two in the phase is due to the grav- 
itational wave frequency being twice the spin frequency 
This step can be performed on the data after down- 
sampling as the rate of change of this phase difference 
will be very low. 

The effect of this extra heterodyne can be seen in 
search for a signal from the Crab in the S2 data. This 
science run of the LIGO interferometers lasted approxi- 
mately two months and overlapped three entries in the 
Crab pulsar ephemeris. The S2 run started on 14 th Feb 
2003, so values of the frequency and spin-down used in 
the initial heterodyning were chosen to be those given in 
the first ephemeris entry prior to the run (15th Jan 2003). 
The second derivative was set to be that taken from the 
ATNF catalogue. The values are shown in Table IIIII 
Once the data were produced the ephemeris values were 



TABLE III: The parameters used in the initial heterodyne 
stage of the Crab pulsar analysis for S2. 




7.31 7.32 7.33 7.34 7.35 



GPS time (sees) 



x 10" 



FIG. 8: The grey points show the phase difference between 
that used in the initial heterodyne and that interpolated from 
a fifth order fit to the ephemeris. The black crosses show just 
the phase difference between the initial heterodyne and the 
individual ephemeris values. 



timing noise removed. There is very little difference be- 
tween the amplitudes for the two cases because the slope 
of the phase difference A(f> are not too steep over the pe- 
riod of S2. However, the extracted value of the phase 



PSRJ0534+2200 



V 


29.8102713888 Hz 


V 


-3.736982 x 10" 10 Hz/s 


V 


1.2426 xl0~ 20 Hz/s 2 


Frequency epoch 


GPS 726624013 



used to calculate the phase given in Eq. (|4.1[) . The dif- 
ference between the initial heterodyne phase and the 5 th 
order phase is shown in Fig. [5] This phase difference is 
used in the extra heterodyne to remove the variation. It 
can be seen in Fig. [5] how a linear fit between ephemeris 
values would be acceptable for these times, with only 
small deviations in phase from the fifth order fit. The 
black crosses in Fig. [5] provide the first step in checking 
the code used for the extra heterodyne stage. The red 
points represent the phase difference used in our extra 
heterodyne step (Eq. (|4.2[) ) to heterodyne each S2 data 
point as calculated using our code, whereas the black 
crosses just show the phase difference between the ini- 
tial heterodyne and the individual Crab pulsar ephemeris 
data points. The fact that these overlap provides a check 
that the heterodyne code is producing the correct phase 
difference. 

If we simulate a signal from the Crab pulsar over the 
period of S2, with parameters /io = 0.5, </>o = 0.0, ip = 0.0 
and l = 7r, we can see how including a timing noise het- 
erodyne step affects the parameter estimation. Fig. [S] 
shows the extracted probability distribution functions 
(pdfs) of ho and (f>o for the signal with and without the 
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FIG. 9: The extracted pdfs for ho and (f>o for a simulated 
signal from the Crab pulsar over the period of S2 with and 
without timing noise removed. 



is strongly affected, mainly due to the the phase offset 
between the start of S2 and the epoch of the initial het- 
erodyne parameters seen in Fig. [5J 

We can simulate a Crab pulsar signal and analyse it 
with and without the timing noise heterodyne step over 
longer periods than just S2 to show its importance. The 
same process as above has been carried out over the pe- 
riod of the S3 run, using the same initial heterodyne pa- 
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rameters and pulsar injection parameters. The extracted 
pdfs arc shown in Fig. [10] and demonstrate that without 
the extra heterodyne the signal is completely lost. The 
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FIG. 11: The extracted pdfs for ho and 4>o f° r a simulated 
signal from the Crab pulsar over the period of S3 with and 
without timing noise removed for initial heterodyne values 
obtained from a fit to the ephemeris over 2003 (see Table WV\ . 



FIG. 10: The extracted pdfs for ho and <f>o for a simulated 
signal from the Crab pulsar over the period of S3 with and 
without performing the extra heterodyne. 



fact that this signal is not seen without the extra hetero- 
dyne may appear at odds with the approximate 2.6 year 
decoherence time stated above, as S3 was only about 8 
months after S2. However, the initial heterodyne values 
used were still those from the Crab ephemeris closest to 
the start of S2 and not those from a more general fit to 
the data over an extended period, as was used to calculate 
Tdocoherence, so it is not the timing noise causing the deco- 
herence in this case, but badly chosen initial heterodyne 
parameters. This does however demonstrate the impor- 
tance of having well-defined heterodyne parameter for all 
pulsars. If we perform a fit to the Crab pulsar ephemeris 
over the whole of 2003 (see Table HVl for fit values), when 
the S2 and S3 runs took place, we can again check the 
impact of the timing noise in S3 (see Fig. ITT]) . We see 



TABLE IV: The parameters of the Crab pulsar for a fit to sec- 
ond order in frequency over the period of 2003 using monthly 
ephemeris data. 
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that, with an extended fit for the Crab parameters, tim- 
ing noise makes little difference over S3, although a slight 
phase offset is present. 



B. Timing noise in other pulsars 

For the majority of pulsars timing noise is most promi- 
nent in the second derivative of frequency, but for mil- 
lisecond pulsars this value is often so small as to be un- 
measurable. For these pulsars, the value of the Ag pa- 
rameter can be used to estimate the cumulative phase 
contribution of timing noise via the empirical relation- 
ship between P and Ag given in Rcf. [32| , Abbott et al. 
fl4| used this as a test of the coherence of the pulsars 
phases over the periods of the LIGO and GEO 600 S3 
and S4 data runs. The generally small values of v for 
the majority of other pulsars suggest that there is not 
significant timing noise. 

There is however, one known interesting target pul- 
sar within our frequency range that has similar proper- 
ties to the Crab pulsar. This pulsar, PSR J0537-6910, 
is young and has the second largest spin-down after the 
Crab pulsar, making it a very interesting candidate for 
our search. It is distant, lying in the Large Magellanic 
Cloud, and has so far only been observed in X-rays. It 
is also a prolific glitcher, showing high levels of timing 
noise [4~0 ] . As dedicated time on the Rossi X-ray Timing 
Explorer satellite is required to time this pulsar, it has 
not been observed as regularly as the Crab pulsar and 
a comprehensive ephemeris does not exist. A regularly 
updated ephemeris would however allow us to track the 
inter-glitch phase and perform an analysis similar to that 
of the Crab pulsar. 



V. CONCLUSIONS 

We have shown how to take account of time delays 
due to the motion of pulsars within binary systems in 
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searches for gravitational waves . This is particularly 
important as a large number of the pulsars within the 
sensitivity band of current interferomctric gravitational 
wave detectors are in binaries. However, as discussed 
in Ref. [l4|, we have demonstrated the importance of 
knowing the binary parameters and their covariances. 

At present our search code contains the three main 
binary models described above. These are sufficient for 
the majority of pulsars, although Tempo contains many 
more models which could be incorporated in the future 
if needed. 

For the cases where large levels of timing noise are 
seen it is important that up-to-date parameters are ob- 
tained allowing us to track the rotational variations. For- 
tunately, the Crab pulsar is constantly monitored in radio 
and therefore has a very well-determined phase evolution. 
Other pulsars in our band are not so well monitored and 
this is especially pertinent if we want to target young 
pulsars (generally the best candidates for gravitational 
waves due to there large spin-down rates), because these 



will generally also be the most affected by timing noise. 
In particular, the young X-ray pulsar PSR J0537-6910 
has high levels of glitch-induced timing noise, and 
an ability to regularly monitor this object would allow 
a more detailed study of this prime gravitational wave 
candidate. 
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